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Abstract 

We investigate the Gubser solution of viscous hydrodynamics at finite density and analytically 
compute the flow harmonics Vn- We explicitly show how Vn and their viscous corrections depend on 
the chemical potential. The difference in Vn between particles and antiparticles is also analytically 
computed and shown to be proportional to various chemical potentials and the viscosity. Excellent 
agreement is obtained between the results and the available experimental data from the SPS, RHIC 
and the LHC. 

PACS numbers: 47.75.-bf, 12.38.Mh, 11.25.Hf 


1 


I. INTRODUCTION 


Relativistic hydrodynamics is a general theoretical framework to describe the collective 
dynamics of high-energy systems near local thermal equilibrium. Its hrst application to 
hadron physics dates back to Landau’s attempt to describe multi-particle production in 

It has become a topic of great interest since the discovery 


hadron-hadron collisions 

of the quark-gluon plasma (QGP) as a n^rly-perfect fluid in the “Little Bangs” at B 


Relativistic Heavy Ion Collider (RHIC) |2|-1^ and CERN Large Hadron collider (LHC) jb-S 
This is supported by the observations that the azimuthal momentum anisotropy of hadronic 


distribution 


10| . characterized by flow harmonics are found to reflect the geometrical 


anisotropy of the overlapping region of two colliding nuclei, and that they are in good 
quantitative agreement with theoretical estimations. Nowadays the viscous hydrodynamic 
modeling is considered as one of the most powerful tools to quantify the QCP medium near 


the crossover phase transition 


ill. 


The recent Beam Energy Scan (BES) experiments at RHIC pose us intriguing challenges 
to study the properties of the medium at hnite density and to explore the QCD phase di¬ 
agram to hnd signs of a critical point 12|. Conserved charges such as net baryon number, 
strangeness and isospin would play important roles in the collisions with lower energies, as 
the differences between particle and antiparticle yields are clearly seen Q, Q- Historically, 
it had long been speculated based on several idealized calculations that the strong coupling 
limit is achieved only at highest energies of RHIC experiments. On the other hand, recent 
improvements in off-equilibrium hydrodynamic modeling motivates us to reexamine the va¬ 
lidity of hydrodynamics in exploring the dense quark matter created at mid-low energies, 
especially since the differential elliptic flow V 2 {pt) is found to remain large in Phase I of the 
BES experiments. The applicability of hydrodynamic models is closely related to the origin 
of fluidity, about which little is known, and thus its verihcation would be a very important 
step towards a full understanding of the hot QCD medium. 

So far many hydrodynamic analyses have been performed numerically because it is gener¬ 
ally quite nontrivial to solve the partial differential equations involved. Analytical solutions 
of relativistic hydrodynamics, on the other hand, can be obtained with certain symmetry 
conditions and they are very instructive in understanding the essence of heavy-ion dynamics. 
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The boost-invariant Bjorken flow 1^ is one such classic example. More recently, Gubser 
found an exact boost-invariant solution of the Navier-Stokes equation which has a nontrivial 


dependence on the transverse coordinate 


Q. 


The latter solution has the advantage that 


one can add azimuthally anisotropic perturbations 181 and a naly tically conipute the 
corresponding flow harmonics Vn including the viscosity effects 18|, ll9| (see, also, 20l-l22|b 


In this study, we investigate Vn at finite density by analytically solving the viscous hydro- 
dynamic equations coupled with conserved currents assuming conformal and boost-invariant 
symmetries. Aside from the fact that the solution itself is new and of theoretical impor¬ 
tance, it gives us a theoretical guidance about the behavior of Vn over a wide range of the 
beam energy for which there are not many numerical simulations 23l-l27l| and the previous 


knowledge obtained through the precision analyses in the RHIC-LHC energy regime, such 
as the value of the shear viscosity r]/s, are no longer fully applicable. We discuss extensively 
the nature of flow in the presence of currents and estimate the beam energy (or chemical 
potential) dependence of The difference in Vn between particles and antiparticles is also 
analytically computed. The results are compared with the experimental data from SPS, 


RHIC and the LHC 


3,[i4, 


29l |. We see that they are in qualitative agreement, which 


suggests that a reasonable description of the low-energy experimental data might be possible 
within a hydrodynamic framework. 

The paper is organized as follows. The basic setup of relativistic hydrodynamics is out¬ 
lined in Section m We then present analytical formulas of the flow harmonics Vn in the 
ideal and viscous cases in Sections m and IIVI respectively. Phenomenological inputs for 
our model are summarized in Section |Vl Using these formulas and input parameters, we 
compare our results with the experimental data in Section IVIl Section VII is devoted to 
summary and conclusions. 


II. HYDRODYNAMIC EQUATIONS 
A. Setup 

We shall consider hydrodynamics of a conformal theory. The system is characterized by 
the local temperature T and a set of local chemical potentials /ij where the subscript i labels 
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various conserved charges of the theory. The flow velocity is denoted by with the nor¬ 
malization u^Ufj_ = —1. The energy-momentum tensor in the Navier-Stokes approximation 
takes the form 

TM- = , (1) 

o o 

where is the shear tensor and r] is the shear viscosity. In ([1]), the conformal equation of 
state e = 3p between the energy density £ and the pressure p has been used. The conserved 
current Jf associated with the chemical potential pi can be written as 

j ^ (2) 

where n* is the charge density and Ki is the charge conductivity. The hydrodynamic equations 
consist of the conservation equations for and 

V^T'^" = 0, V^Jf = 0, (3) 


where is the covariant derivative. 

Since there is no intrinsic mass scale in a conformal theory, the energy density e and the 
charge densities n* can be generically written as 


e = T^f 


VT’T’"' 


TLi piT Pi 


(t!l t!l 

VT’T’"' 


( 4 ) 


With a view to applying to heavy-ion collisions, we shall focus on the following representative 
situation. We assume that there is the leading current = nu^ + ■ ■ ■ (‘baryon number 
current’) and the corresponding chemical potential p is treated to all orders. In addition, 
there is one subleading current = hu^ + ■ ■ ■ (‘isospin number current’) whose chemical 
potential p is small and treated only to linear order. We take p to be ‘orthogonal’ to p, in 
that e{p,p) is invariant under a sign flip p -H- —p (i.e., cross terms like ppT‘^ are absent). 
With these assumptions, we can parameterize 


e = T^f , n = pT^g , h = pT^g . (5) 

The last equation may be written as h = px where x oc susceptibility. 
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B. Gubser flow 


We shall solve the hydrodynamic equations 

2tx I 


= cosh 


tanh 


-1 


2^2 _|_ ^2 2;2 


for a given flow velocity 
2txj 


= sinh 


tanh 


-1 




( 6 ) 


and = 0. The parameter L is the characteristic length scale of the system. In 

heavy-ion collisions, it is roughly the transverse size of the colliding nuclei. Eq. (|6]) is called 
Gubser flow Q, 12 1 expressed in the coordinate system 

ds^ = —dr"^ + r'^dC'^ + dx\ + x^dcfP' , (7) 

where r = — x| is the proper time, C, = tanh“^ ^ is the spacetime rapidity and x± = 

xl + X 2 is the transverse coordinate. The condition = 0 means that the flow is boost 
invariant along the beam (xa) direction. 

Gubser flow takes a very simple form in a cleverly chosen coordinate system x^ which is 
related to the Minkowski coordinates via a Weyl rescaling of the metric. 

ds^ 


ds = 




= —dp^ + cosh p{dQ^ + sin^ Qdcj)'^) + d^^ , 


where 


sinh p = 


+ x] 

2L7 


tan© = 


2Lx I 


IP ^2 _ 2;2 


( 8 ) 


(9) 


In this coordinate system, the flow velocity is simply it^ = S^. In addition to the boost 
invariance, the flow respects the 0(3) symmetry with respect to the ‘polar’ angles (0,^). 
Variables in this coordinate system will be denoted with a ‘hat’, e.g., it^, e. 

III. INVISCID CASE 

In this section, we solve the hydrodynamic equations (E]) in the ideal case p = Hi = 0. 
We then deform the solution in the azimuthal direction 0 and compute flow harmonics Vn- 


A. Isotropic ideal solution 

The isotropic solution (i.e., independent of 0) has been obtained already in Q, Q in 
the presence of a current = nu^. Assuming that all the quantities depend only on p, we 
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can readily solve the hydrodynamic eqnations for eo and ho in the coordinates ([8]). We then 
perform the Weyl transformation back to the Minkowski space Eq = Eq/t^, uq = uq/t^ to 
get 


^0 = ToV 



1 

r"^(coshp)®/3 ’ 


( 10 ) 


m = 1‘oT^g ( K ^ ; 

Jo/ cosh p 


ho = poT'^g ( oc ^ 2 
Jo/ T-^ cosh p 


( 11 ) 


These equations can be solved for To and Pq. It is consistent to look for the solution where 
To and Pq have the same p-dependence such that the ratios a = Po/Tq, a = Pq/Tq are 
independent of p. We find 

C aC 


To = 

and therefore, 


^0 — 


r(coshp)T3 ’ 

/(o)C‘ 


ho — 


r‘^(coshp)8/3 ’ r^cosh^p’ r^cosh^p' 

The parameter C is related to the particle multiplicity to be extracted from the experimental 
data. For a massless particle species i (‘pion’), the relation is jlbl. Il9| 


r(cosh p)T3 
ag{a)C^ 


Po — 


no = 


aC 

r(cosh p)T3 
q;p(ck)C^ 


( 12 ) 


(13) 


dNi 


9i- 




dY TT 

where Y is the momentum rapidity and pj is the degeneracy factor. 


(14) 


B. Anisotropic ideal solution 


We now perturb the solution anisotropically to introduce the cos n0 dependence. In doing 


so, we shall focus on the early time regime r <C T (or p —)■ —cx), see ([2])). As observed in jl9| . 


in this regime the perturbed solution is fully under analytical control including the viscous 
case to be discussed in the next section. 

Following 17|], we consider the following deformation of the isotropic solution 


Eo ^ E — do(l — EnA-dY , 

no ^ n = ho(l - EnAd'f , 

u^) = (1,0,0,0) (l, -EniEsg^^deA, -eni'sg'^'^d^A, O) , (15) 
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where 




2Lx± 


X, 


cos n 0 , 


( 16 ) 


is proportional to the spherical harmonics 0 ) + 0 ) in the early time regime 

T L. Note that we preserve boost invariance = 0 in this paper, but the case n'’ 7 ^ 0 


was also considered in 


17| . e„ is the eccentricitjll] which we assume to be small e„ 1 


and keep only linear terms in e„. 6{p),6'{p) and Usip) have to be determined by solving the 
hydrodynamic equations linearized around the isotropic solution. Plugging flThl) into ([3]), we 
hnd the following equation for 6' 

(18) 


dpS' = 


- -^n{n + 1 ) 

3 cosh^ p 


This turns out to be exactly the same as the equation satished by 5 


17l |. Therefore, in 


the ideal case we have 5 = 5', which means that Tq and Pq are rescaled by the same factor 
T = To(l — enAS), p = po{l — enA6) and p = po{l — enA6). The ratios p/T = Pq/Tq = a 
and pjT = po/T^ = a are thus unchanged. At early times p —)■ — 00 , the right hand side of 
firsj) is negligible and we can set 5 = 1 [19 1. 


C. Vn Sit finite p 


In order to compute flow harmonics we use the Cooper-Frye formula 


30| 


(2vr) 


dNi 


{-p^da^) ( exp 


u ■ p + kpt 


+ (5/ oc 1 + 2vn{pT) cosn(j)p ,(19) 


dYpxdpTdcjyp V V ^ 

where we assumed the Boltzmann distribution and 6f is the deviation from the equilibrium 
distribution. The use of the Boltzmann distribution (rather than the Fermi/Bose distribu¬ 
tions) may be justihed for the purpose of computing the integrated [l9j. pt generically 
represents a set of chemical potentials for net baryon number, isospin and strangeness. We 
assign k = ±1 for particles with positive/negative quantum numbers mentioned above, and 
k = 0 for neutral particles with respect to the corresponding quantum number. In princi¬ 
ple, since we are assuming conformal symmetry, the formula (IT^ should be used only for 

^ In a conformal theory, the definition of eccentricity requires some care. We use 17 [Q 

J cos n(j) 




( 17 ) 
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massless particles, or particles that can be approximately treated as massless (i.e., pions). 
However, for the sake of discussion in Section IVl Cl we shall later introduce massive particles 
and compute their Vn in the ‘probe approximation’, namely, by neglecting their backreaction 
to the flow velocity. Since we add in particles in the final state that do not exist in the fluid, 
the total energy is not conserved at freezeout. But the fraction of the change 6e/e ^ q-^It 
is exponentially suppressed by the mass m and will be neglected. 

The integral in flT^ is taken over the hypersurface E of constant energy density where 
the kinetic freezeout occurs. In the ideal case, constant e means constant T since a = fi/T 
is a constant. Let us write the condition of constant energy density as 


^4 o4 

£(r, IX. 4>) = T‘f{a) = ■ 


( 20 ) 


Typically, Sc is of the order of the critical energy density of the QCD phase transition. We 
take Ec = 1 GeV/fm^ in this paper. Following 19| , we assume that the condition fl20|) is 
reached within the early time regime r L where we can use the approximate solution 
(1T5|) . The parameter B in (l20ll is then related to the (position-dependent) freezeout time r/ 
as 

(2L)5 ^ f 2Lx I " 


Tf{x±,(j)) = 


H3(L2+xi)2 




L2+x2 


cosn0 


( 21 ) 


For consistency with our early freezeout scenario, we must have 5^ S> 1. 

Under these assumptions, the integral flT^ can be performed analytically and the inte¬ 
grated Vn is obtained from the formula 

/‘‘PrJi, ■ ' ' 

In the ideal case 5f = 0, rv does not depend on k since the factor = e*’" cancels in the 


ratio 


6 rj 


The result is 
9 F(3n 


19 


64 F(4n) 


128\ ”■ 2 ^^(3ri + 2)^(n — 1) 

^ V2/ 2(4n + 1) 


B 


—3n 


(X 


/3U dN 


(23) 


This determines the a = p/T dependence of Quite generally, f{a) is an increasing 
function a. On the other hand, dN/dY is a decreasing function of a. We shall see that, in 


2 See (61) of 19|. We have corrected a mistake by a factor of 2 in the overall normalization. The same 
comment applies to (l48l) below. 





















heavy-ion collisions, the latter dependence is stronger, and as a result (|23|) is a decreasing 
function of a, or equivalently, an increasing function of the collision energy y/s. Incidentally, 
we note that the directed flow Vn=i vanishes, consistently with our assumption of boost- 
invariance. 


IV. VISCOUS CASE 

We now turn to the viscous case 77 , Ki 7 ^ 0. Although the system is out of equilibrium, 
from the Landau matching condition we can dehne the local T and p using the same relations 
as in equilibrium 

e = T^f , n = nT‘^g , n = flT^g , (24) 

but now /i/T cannot be a constant. 


A. Isotropic viscous solution 


First consider the isotropic case = 5^. Although p/T in ([2]) is not a constant anymore, 
it depends only on p (see below). Then we still have = n6^ so that 

1 


n oc 


is the same as in the ideal case 


cosh^ p ’ 


(25) 


17|. However, the solution of the Navier-Stokes equation ens 


has an extra p-dependence proportional to the shear viscosity ru In the case of vanishing 
chemical potentials, this p-dependence can be obtained exactly |l 6 | 


^NS = 


fC^ 


(coshp)®/^ [ 


1 + 


sinh^ p 2 F 1 


3 7 5 
2 ’ 6’2 


; — sinh^ p 


n 4 


(26) 


where fj = is independent of p. 

However, at hnite density, fj will depend on p, and this makes it difficult to hnd an exact 
solution. Related to this, p can now depend on both e and n, and this relation can be model- 
dependent. We can get around this problem by assuming that p is small. Specihcally, we 
rescale p by the entropy density s 


Tj = Tj 



s, 


9 


(27) 









as is often done in hydrodynamic simulations. We then regard 77 as a small parameter 
iv rNJ 0(10 ^)) and keep only terms linear in rj. In this approximation, we may replace /i/T 
and s in fl27|) by their equilibrium values at 77 = 0, namely, fi/T = a and 

i 

-h{a). 


1 

-{eo+po- pono) = — - 

j-o T'* cosh p 


4 

-/(a) - a^g{a) ) = 




We then find the solution valid to 0{ri) 


ens = T‘f 


/(tt)C- 

T'*(cosh 

/(o)C* 


uns = pT^g (^) 


r'^(cosh 

ag{a)C^ 


^ Ah(a)p(a) . , o ^ 
2h{a)fj{a) ( 


cosh^ p 
3 7 5 


(28) 


2 ’ 6 ’ 2 


; — sinh p 


1 - 


f{a)C 


2 7 


(29) 


(30) 


cosh^ p ’ 

where in the second line of fl29p we focus on the early-time regime where p is negative and 


larger 

Using 


and 


, we can eliminate C 


[a 




fiP 


f(a) 1 - 


2h{a)T]{a) f 

/WC 


m 


Writing 


^ = « + ^«(P) , 


we find the deviation from constancy due to the viscosity 


6a{p) 


2/777 


3g 


hp /' L"^ + x]_ \ 

77 J 


where 


7(a) = 


A _L M 
3a 3g / 


Note that 7 ( 0 ) oc a as a —)■ 0. At the freezeout time r = tj, we have the relation 

y Cr A freezeout ~ O A pK ^ 2 A)^ ' 

Finally, we can solve for T and p using ((321). The result is 


T = 


C 


r(cosh p)2/3 


a 


p = 


aC 


r(cosh p)2/3 


1 + A| -_iL 

a 


(31) 


(32) 


(33) 


(34) 


(36) 


(36) 


The viscous Gubser solution is known to become unphysical (the temperature becomes negative) as 


—00 


16l |. Physically, this corresponds to very early times and/or very large values of a;_L. Our results 
are not sensitive to these regions. We can simply choose the initial time of the evolution to be small, but 
not too small. Besides, all the a:j_-integrals to be^Performed below are fully convergent at a:_L —>• 00 . 




























B. Anisotropic viscous solution 


We now perturb the solution as in (fT5|) . First consider the current in ([2]). p/T now 
depends not only on p, but also on 0 and 0. However, the dependence is of order p. (See 
fl43p below. Remember that for the ideal solution p/T is constant even in the anisotropic 
case.) Therefore, if we neglect terms of order 0{Krien), we can approximate ~ nu^. Then 
flT8|) is still valid and we get 


n = fiT g[- 


fi\ ag{a)C^ 


cosh^ p 


(1 - enAf 


(37) 


As for the energy density, we hnd 


£ = 



where 



f{a)C* 
r^(cosh p)^/3 


/ _ 2h{a)T]{a) 

\ f{cy)C 




5 


1 + 


h{a)i]{a) 

2f{a)C 



(38) 


(39) 


From the constant energy condition 


C^B\ 

£ = T:rvTT/(a) = £c, 


(2L)^ 


we can determine the freezeout surface T{x±,(j)) in the viscous case 19| 

{2Lf ( 2Lxi 


r/(xx,0) ■53(^2 + 3.2)21^- 2 ( 2 L )4 "^"\ L ‘^ + xI 

where the ‘Knudsen number’ is proportional to the shear viscosity 


(40) 


COS ncf) I , (41) 


K = 


f{a)C 


and 


are also modihed as + = « + 5a' where 


5a — 5a{\ T — 5a | 1 T C: 


2Lx I 




cos ncj) I , 


and 


T = 


h = 


C 


r(cosh p)A3 
aC 

r(cosh p)A3 


^ , 5a f 1 g' 

a g 


5a f 2 g' 
l-enA+ — { - -- 
a g 


(42) 


(43) 


(44) 
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C. Vn Sit finite // and rj 


The computation of Vn is more complicated than the /r = 0 case. This is because e = const 
does not mean T = const, and therefore one cannot treat T in the Boltzmann factor (IT^ 
as a constant when integrating over the hypersurface of constant energy. In order to cope 
with this, we write fl44p as 


T f 

T rri OJ r 

= Tc- 


(45) 


where 


Tr = 


CB 


— Tn 1 — CnA — 


hr] f + Xj 


2 / 3 ^ 


2L M " 2/C V 2Lr 

is constant by virtue of (140]) . We then expand the Boltzmann factor ac 


(46) 


exp 


u ■ p + kjj, 
T 




, u ■ pTof . , 


|i + r + (1 + £„.4)| , (47) 


where we approximated ~ To(l — CnA) in the 0{6a) term. 

The first term in fl47l) . proportional to unity, gives the same result as in ll9P 
9 r(3n) f 128\"' ^2 f'^\ f n'^ {3n + 2)'^ (n — 1) 3n^{n — l)K 


64r(4n) \B^ 


2(4n + 1) 


16(3n- 1) 


(3n^ + 3u + 2) K48) 


Note that = 1 — 0{nK) for u 3> 1 (see, however, 3^). The second term in fl47)) 

leads to a new order 0{K) contribution to To compute it, we borrow some results 


from 


19|. First, the perturbed flow velocity on the freezeout surface has the following 


components in the coordinates ([^ 


ui_ = Mo± + cos ncj), u^ = Su^Cn sin ncf ), 


(49) 


^ In this subsection we set /i = 0. The case with fi / 0 will be treated in the next subsection. 

® For simplicity, here we ignore the contribution from the nonequilibrium part <5/ in (1191) . This has been 

computed in Q for a particular choice of Sf. However, its n-dependence is strongly affected by the 

choice of 6f which is not unique 31]. Moreover, even the overall sign of this contribution is sensitive to 

the pT-cutoff. 
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where 


u±o = 2x± 


(2L)- 


2 13 ’ 


5u I = 


= 


2Lx± 


B^{L^ + xl) 

3(2L)5 

B^{L^ + xiy \L^ + x 

3n {2Lf 


n—l 


2 
± 

2Lx± 


L{n{L‘^ — x^) — ^x\_) , 


( 50 ) 


* 2 + xlY \V‘+ x] 

(The viscosity can be neglected here.) The exponential factor in the Boltzmann distribution 


reads 


p ■ u 1 


—rriT cosh(C — F) + PtU± cos(0 — 0 p) — sin(0 — </>„) 

x± 


= U + eJ[/,(51) 


where = \/is the transverse mass. The volume element of the constant energy 
hypersurface is 

( ^ 7 _ ^' 7 “ \ 

rriT cosh(C — V) — Pt cos(d) — (i>„) 7-^ H —— sin(0 — (bp) — -f ) dCdx^dd) 

OX± X± 0 (j) J 

(52) 


where Tf is given by fl4T|) with the viscous term set to zero. Finally, we need the more precise 
version of (l35l) 

(53) 

Armed with these formulas, let us decompose the contribution from the second term in 
fH71) as 

= ((5Ji + (5J2 + 5J3)enCosn0p, (54) 


corresponding to the three terms in (1521) . Consider 6Ji first. To 0{en) we have to evaluate 
6Ji 


~ Cn—^ 77 ^- rriT [ dx±x± [ d(d(j) cosh.{(— Y)e^6U (+-^ + k \ . (55) 


53 


4/ 4/ 


This can be efficiently evaluated using the trick introduced in 19| (see Eq. (73) there). The 
0-integral gives Bessel functions In{z) where 

Ptu±o 2xxPr(2T)3 


2 ; = 


Tr 


TpB^{L^ + xl) 


2 'iS • 


(56) 
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This can be expanded as In{z) ~ z" anticipating that the subsequent p^-integral is domi¬ 
nated by the region z < 1. We thus find 


, ^ 2L-fKe^^ , Attz 

dJi ^ - — - rriT / dx±x± 


X 


S3 

L 

4 / 


2"(n - 1 )!m±o 
TUt 


6u 


x± 


nK,{mT/T,) - ^ (i^o(mr/T,) + K,{mT/T,))) + kKrimr/T,) 


2L'yKe’^°‘^ ^2 /Gfip-r \ " n(n — l)r(3n) 

S3 ^ \T^B^) (3n - l)r(4n) 

T (^nKiimrlK) - (K„(mT/X) + K^XfiiTlX))) + kKi(mT/T,) 


■ (57) 


The correction to can be calculated from the formula (cf. fl22|) l 


r 1 In dprPT^Ji er, 

ov„ = 


dprPrJ^ 2 

where is the azimuthally symmetric part (cf. Eq. (45) of jl9j |) 

16L3 


(58) 


= ATirriTKiiniT/T()- 


^kOL 


S3 


(59) 


In the massless case ttit = Pr-, the integral can be done exactly and we find 

dvj ^ 97^ / _ ^ n^(n- l)r(3n) fn \ /n\ 

en 128 V Af)\B^) (3n-l)r(4n) V2 / \2J 

and from fl2^ . 


(60) 


fo). pTf 3/7 ri(ii +2)(4ii+1) 

yideai 4 \ 4// (3n — l)(3n-|-2)2 

It is important to emphasize that fl6T|) is induced by the combined effect of the chemical 
potential and the viscosity. It vanishes when 7 = 0 or a = p/T = 0 because 7(0) = 0 
(cf. fl3T)) b Compared with fl48|) which schematically reads ~ —nK, we notice 

that fIM]) is not enhanced by a factor of n, hence subleading at large n. However, it is 
the leading contribution to the difference in Vn between particles {k = 1) and antiparticles 
(/c = —1). If/i = /iB>0is the baryon chemical potential, the protons have larger than 
the antiprotons. We shall study this effect in detail later. 

In fact, for protons the approximation niT ~ Pt is not valid. Instead, we now assume 
mT S> T and reevaluate Note that when tut ^ T, 5Ji is parametrically larger than 
(5J2,3, so it is enough to consider only 5Ji. 
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hT, 


When rriT S> T, the Bessel function is independent of the order 

Ki{mT/Tc) « 

so that flFr|) becomes 


2mi 


iLg-my/Tc 


27 /^ 6 ^“ 3 (Upt\ n(n-l)r(3n) 


On the other hand, from Eq. (47) of 19j, 


^(n-^\+k 




/ 64^ Y 


n 


1 ) 


2(3n + 2) 


(62) 


(63) 


(64) 


477, + 1 

The pT-integral can be evaluated by the saddle point at = yJnmTc for m 3> nT^ and we 
obtain 

dv^ yit' 77 ,( 477 , + 1) 


ddeal 


I Ti TTl 

4 (3ij - l)(3n. + 2) V47V2~fr'^*' 


(65) 


The fc-independent part is order ^ nK, but we shall see later that it is numerically 
small for realistic values of m because the factor jf/f is small. The /c-dependent term is 
again of order 0{K) without an enhancement by a factor of n. 


The evaluation of (5J2,3 in llhTl) can be done similarly, though it is considerably more 
tedious. Here we only show the final result in the massless case 777 = 0, relegating the details 
to Appendix 


97/7 


k- 


3f 


128 

'W 


77^(77 — l)r(377) fn 


(3n- l)r(4u)'^ I2 ^ V2 


77 


( 66 ) 


en 128 V“ 4/ 

Comparing with (legD, we notice that 5v\ = Actually, this relation was repeat¬ 

edly observed in 19| when computing other contributions to We do not have a simple 
explanation for this. 


Summing all the contributions including the previously computed term 19|], our final 

result of the viscous correction 5vn in the massless case is 
5vr,. K r(377) f 128\ " /77\ 77^(77 — 1) 


K r(3u) /128y , 
256r(4u)V53^ V2 


3u - 1 


27 


X • 


— (377 377 2) 97 — -M 


377 



(67) 


4 ' V 2 

The second term in the curly brackets is the new contribution at finite density. It is sub¬ 
leading in 77 , and actually the factor 7 (a) is also numerically small. However, it gives the 
leading contribution to the difference in between particles and antiparticles. 
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D. Isospin chemical potential 


In the previous subsection, we computed of particles which couple to the ‘large’ chem¬ 
ical potential fi. Here let us compute Vn of particles neutral under p but charged under jl. 
We have in mind the charged pious in the presence of the isospin chemical potential. We 
start with the formula (cf. ([37j)) 


_ Oig{o!)C^ , .3 ^ _ 2 ~/ /rji\ 

r2 gyn/T). 

cosh p 


( 68 ) 


We treat d = p/T as a small parameter and keep only terms linear in a. Dividing by 
from (Hill and using p/T = a + 6a{l -|- e„^), we hnd 

ag{a){l — SnA)^ _ /I 


h 

T 




a \ 1 + 5a 


- + ^ 

a g 


9 


(1 -|- SnA) ) . (69) 


The fugacity factor thus becomes 


gtp/T S3 e‘“ ( 1 + kaSa ( ^ ^ J ) (^ + ^ 


(70) 


As before, the factor e^" drops out in the computation of Vn- We see that the only difference 
from the previous case (l47D is that k is replaced by 


a 


a 


1 + 


ag' ag' 


9 9 

Thus the hnal result is the same as fl67)) except that k is replaced by fl7T]) . 


(71) 


V. PHENOMENOLOGICAL INPUTS 

This section serves as a preparation for the next section where we compare our results 
with the experimental data. 


A. Models 

In order to make quantitative predictions, we need models for the functions /, g, g dehned 
in (|5]). Here we consider two extreme scenarios in terms of the interaction strength. 
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1. Free quark-gluon gas 


The energy density of free, massless three flavor QCD is 


Uq = 


dp 

Opq 


— PqT I 1 + 


. = 3p=^T‘ + 6 5: 

q=u,d,s 

f^q 

7r2T2 






T^ + 


+ 


Svr^ 


(72) 


where Pu = ^ + ^, ^ ^ and /i^ = ^ — Ps- Fb, I^i and ps are the baryon, 

isospin and strangeness chemical potentials, respectively. Since the net strangeness is zero 
in heavy-ion collisions, we set /is = /Ws/S and obtain (a = Pb/T) 


19 


a 


f(a) = —TT^ H-h 

^ 12 3 


a 


Sdvr^ ’ 


(73) 


9(a) = 9B(a) = 2 ll + , g(a) = gi(a) = ^ ■ 0^^) 

It turns out that, due to the large denominators Qvr^ or Svr^, the effect of g and g on is 
numerically small. 


2. Af = 4 SYM at finite TZ-charge chemical potential 


Next we consider strongly coupled A/" = 4 supersymmetric Yang-Mills theory at finite 
i?-charge chemical potential /i. This theory is conformal, and in the limit of strong cou¬ 
pling and at large Yc, the thermodynamic quantities can be computed from the AdS/CFT 


correspondence. The results are 


33| 


e = 3p 


8 ¥ 



2/i2 

37r2T2 



3 



3a 1 + 


2/i2 

37r2T2 



(75) 
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2/i2 

37r2T2 



2 


(76) 


where n is the i?-charge density. The shear viscosity is given hy p = ^. 

There are uncertainties when treating this model as a proxy of strongly coupled QCD, 
such as the value of and the proportionality constant between p and pb- However, in 
practical fits, the normalization of / oc can be absorbed by a change in L (cf. Eq. fl20l) b 
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Moreover, as long as /i ~ 0(ij,b), the two functions fl73|) and d75|) are qualitatively not so 
different in shape for fi ~ 0{T). As a result, the quality of hts is similar in the two cases 
despite the huge differences in the underlying dynamics. Therefore, in the next section we 
show only the results based on fl75D and flTlD . 


B. Preezeout conditions 


We employ the following phenomenological parametrization 


34j | of the freezeout temper¬ 


ature T and chemical potential hb (in units of GeV) as a function of the collision energy 
y/s (per nucleon, in units of GeV) 

d 


T{hb) = a — bjj,^ 




iB ) hs — 1 , /- ) 

1 + ey/s 

with a = 0.166, b = 0.139, c = 0.053, d = 1.308, e = 0.273. This gives /is/T as a function 
of ^/s as shown in Fig. [H The curve is well approximated by /i^/T d/{aey/s) ~ 29/i/s. 
Actually, T and /is here are the chemical freezeout parameters which are in general different 
from those entering the Gooper-Frye formula (IT^ used at the kinetic freezeout. However, in 
our model only the ratio /is/T matters, and this ratio is roughly constant as we have seen. 
We thus use the relation in Fig. [T]for the evaluation of 

The parameter C also depends on_^ via flTT)) . We use the following empirical formula 
for the charged particle multiplicity 351 

dN^h 


(77) 


4C^ 

2_ =3 148(yi)“ 


(78) 


dY 71 

where the factor of 2 counts the degeneracy between 7r+ and vr”. From (HOI) and fITSll . we 
see that the Knudsen number fl42D behaves as 
_ ShfjBViF, fjB 

~ 377rs/2(^)o.s ~ /v2(yi)o.3 ~ ^nrw^Fff’ 
where L is in units of fermi and = 1 is in units of GeV/fm^. Putting aside the potential 
/is-dependence of r) = r^/s, we see that K is an increasing function of /is (up to /is vrT 
in our model) or a decreasing function of i/i. 

In fact, up to the RHIG energy, we hnd that the following parametrization also gives a 
good description of the data [39 1 

dNch ^ ^2 In . (80) 


dY 


1.41 


18 















FIG. 1. /is/T at freezeout as a function of the collision energy ^/s (per nucleon). 
We shall also use this in Section IVIBl 


VI. COMPARISON WITH THE EXPERIMENTAL DATA 

In this section, we compare our results with three different experimental data: (i) the n- 
dependence of measured at the LHC; (ii) the collision energy dependence of V 2 measured 
at the SPS; (iii) the difference in V 2 between particles and antiparticles measured at RHIC. 


A. Higher harmonics Vn 


The CMS collaboration at the LHC has measured the Pr- integrated in lead-lead colli¬ 
sions at y/s = 2.76 TeV up to rather high orders (n < 7) [3^. Using fl2^ and fl^ together 
with the phenomenological inputs in the previous section, we can evaluate and compare 


with the CMS dataE The result is shown in Fig. |2l Here we set e„ = 0.018 for all dif¬ 


ferent values of n. Taking e„ to be independent of n may be a good approximation for 


® The CMS uses the pr cuts 0.3 < pr < 3 GeV while our analytical result is integrated over all pr- We 
checked that the quality of the fit is unchanged by introducing cuts in our model. 
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FIG. 2 . Vn versus different values of n from 2 to 7 measured by the CMS collaboration ( 0 - 0 . 2 % 


centrality) 


36l |. The black solid curve represents the ideal hydrodynamic result, while the blue 


dotted curve and red dashed curve correspond to the viscous results with rjjs = 0.08 and r]/s = 0 . 2 , 


respectively. See, also. Ref. 


3 ^. 


the very central (0-0.2% centrality) nucleus collisionsl3 The parameter L is set to 17 fm. 
The corresponding value of B in f|48ll is ^ 26.7 which is consistent with the assumption 
B^ > 1 . 

As a matter of fact, since ~ 0 at the LHC, the new term at fiB > 0 (the term 
oroportional to 7 in fl67)) ) is negligibly small, and the present £t could have been done in 
19l | treating i? as a htting parameter. By expressing B in terms of observables as we have 
done here, we can test our result at lower energies or higher chemi cal p otentials fis ~ 0{T). 
Note that since B^ is larger at lower energies, Vn ~ jrj) j^gj fjeQj-gaggg faster with 

n, and this will make the measurement of higher harmonics difficult at low energies 


B. Energy dependence of 


Next we turn to the energy dependence of the elliptic flow Vn =2 for which there are already 
a wealth of experimental data from the SPS and the RHIC BES program 28|, l38|]. We 


^ The value 0.018 may seem a bit too small. This may be due to our nonstandard definition of (HZD. 
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compare our formulas (H 8 |) and (l67|) (with A; = 0) for n = 2 with the SPS, mid-central data 


collected in the low energy region ^/s < 20GeV [281, The result with three different 

values of r]/s is shown in Fig. [3]where we tried both fl78l) and flSOll . the latter actually gives a 
better description of dNch/dY in this low energy region. The other parameters are chosen as 
L = 15.5 fm and 62 = 0.32. The value of L here is slightly smaller than the one (L = 17fm) 
used in Fig. [2l This is consistent with the perception that the QGP droplet is larger at higher 
energies at the time of thermalization. The rise of V 2 with energy is nicely reproduced by 
our formula and attributed to the rise of dNch/dY. It turns out that the newly calculated 
viscous correction in Section HV Gl (the last term in (1671) ) is numerically very small (about an 
order of magnitude smaller than the hrst term in (1671) ) even in the highest density region. 


Unfortnnately, this £t, which agrees reasonably we 


the high energy RHIG data at y/s = 200 GeV 


40 


1 with the low energy data, overshoots 
4l| in similar centrality bins by a factor 


of 2 (assuming that 62 is independent of energy). This is because the rise of dN/dY with 
energy is too steep. If we artificially reduce the exponent i n (j^ as 0.3 —?■ 0.23, for example, 
we get a decent description of V 2 over a broader range in \/io Alternatively, the dependence 
V 2 ~ {dN/dYY from (123|) may be too strong, and the experimental data actually suggest 


a weaker diV/dF-dependence 


While we do not have a resolution of this problem in 


the present framework, it seems qnalitatively correct that V 2 is directly proportional to the 
multiplicity to some positive power, and therefore it is an increasing function of -y/i (see, 
also. Section VII of 19|). 


C. Difference in Vn between particles and antiparticles 


Finally, we investigate the difference in V 2 between particles and antiparticles which has 


been measured by the STAR col 


from theoretical viewpoints 


42 


aboration at RHIG [13 


, and attracted some attention 
45| . For a hadron with the quantum numbers 


(baryon number, isospin, strangeness), we assign the fugacity factor 

B/jb 1 I SjJs 

exp ■ 


(81) 


® The SPS data do not have a low-pr cutoff while the RHIG data have pt > 0.2 GeV. Our analytical 

formula, integrated over all pT, should fare better with the SPS results. 

® We thank Anton Andronic for correspondence about the SPS data. 

Note that (l78l) is for central collisions. The exponent may indeed be smaller for mid-central collisions. 
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FIG. 3. The energy dependence of V 2 compared with the SPS data 
We used (1781) and (fHOjl in the left and right plots, respectively. 


for mid-central collisions. 


(S' = — 1 for the strange quark.) In heavy-ion collisions, fij < 0 since the colliding nuclei 
are neutron rich, and fis ~ hs/S since the net strangeness vanishes. The latter condition 
implies that we should not treat u.g as a small perturbation. Indeed, various estimates of 


rs = based on the SPS and RHIC data, and also from lattice QCD 


all found similar values within the range 0.21 < rs < 0.27. We thus regard fis as a shift of 
fiB for strange hadrons and treat it as a htting parameter, anticipating that the value of r 5 
should come out in the window 0.2 < rs < 1/3. On the other hand, we regard r/ = fij/fiB 
as a small parameter compared to unity and use the result obtained in Section IIVDI 
Let us dehne the difference in between hadrons X and antihadrons X as 




(82) 


This can be evaluated from flH7|) and fl7T]) . Focusing now on the elliptic flow case n = 2, we 
can immediately write down the following ‘master formula’ 


. V 6144 

Xvo = 62 - rlK 

^ 35R6 ^ 


e + rsS + r,/|l + ^-“®' 


(83) 


9 9i 

By construction, (l83|) has been derived for massless particles. In the massive case, we observe 
that the following ratio 


X,ideal 

V2 


40 


B + rsS + riI[l + 


ag agj 


9 


9i 


(84) 
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is exactly independent of mo This is due to the nontrivial cancelation of pr-integrals such 
as fl58l) in the ratio for the /c-dependent part. In order to get Au 2 itself, we must multiply 
dBl by [igl 

^ , 2“ J dprpj. (^K,(mr/T) + 4K„(mr/T)} 

^ ^35B^T f dpTPT'^TKi{mT/T) 

The m-dependence of flH^ is sensitive to the cutoffs of the pT-integral, but overall the 

dependence is not very strong. For simplicity, in this study we ignore the m-dependence of 
^ideai flH^ for all hadron species. It is not difficult to implement this mass effect, but 

there are other subtleties which are not taken into account, either o Clearly, it is desirable 
that the experimental results are plotted in the form (IMj) in order to avoid various systematic 
uncertainties. 

The most important feature of (l8^ or fl84|) is that is proportional to both the 
shear viscosity r] and the chemical potentials. (Remember that 7 as dehned in fIMl) is 
roughly proportional to fiB-) This in particular means that AnJ can be nonzero in viscous 
hydrodynamics in the presence of the isospin chemical potential. 

Let us confront fl 8 ^ with the data. The STAR collaboration has measured Av^ for 


14| . This is plotted in Fig. 0] together with our ht based on 


A = 7r+,K+,p,A,E- [13 
with rj/s = 0.2. We have used L = 15.5 fm and 62 = 0.32 as in Fig. [3l and used the ht pa¬ 
rameters rs = 0.23 and r/ = —0.15, the former is consistent with our expectation mentioned 
above. The steep rise of Av^ for baryons towards the low-i/s region is due to the rough 
proportionality Av 2 oc 7 oc ps. Compared to this, the dependence of the factor l/R® 
is subleading. Since p. A, JF+ have = (1,|,0), (1,0,—1), (1,—|,—2), (0, |,1), 

respectively, we expect the ordering Auf > Av^ > Av^ > Av^^ > 0 for reasonable val¬ 
ues of > 0 and rj < 0. This tendency is obeyed by most data points except a few in 
the low energy region. We note that the data point at ^/s = 11.5 GeV should not be 


taken seriously because, according to the STAR collaboration [^, this data point is afflicted 

with ‘additional systematic effects which are not included in the error bars’. In Fig. IH we 

As already noted in Section fill Ci we introduce massive particles in the probe approximation, namely, we 

let these particles flow with the same flow velocity and neglect their backreaction to the velocity. While 

this causes some inconsistencies such as energy nonconservation, we expect that the essential features of 

(|84l) (the proportionality to the viscosity, /i’s and the corresponding quantum numbers) are robust. 

For instance, the STAR collaboration uses 0-80% centrality events to measure Av 2 - This reduces the 


effective value of by about 20% 
computing Av 2 - 


14| and i2)3tly cancels the above mass effect for baryons when 














have also included our prediction for the ff-baryon. Since has S = —3, we expect that 
Av 2 ~ /x_b — 3/i5 is smaller than other baryonsl^ 

Concerning the pions, the negative can be naturally explained by the negative 

isospin chemical potential. However, the magnitude is problematic. Our choice rj = —0.15, 
which describes the pion data very well, is too large compared with the value rj —0.02 ~ 
—0.03 extracted from the SPS data j^, 47\. We may dial rj down to, say, rj ^ —0.1 without 
spoiling much the quality of the £t, but not further down. On the other hand, the 

other hadrons (p, A,S“, JW) are more or less unaffected by rj and can be well htted eyen 
with rj = —0.02 and rs ~ 0.2. This may be an indication that there are other mechanisms 
to generate the difference Av 2 which predominantly act on the pionsUl 

In the large--y/i region, our result tends to slightly overestimate Av^■ This is partly due 
to the too fast rise of with energy as mentioned before. However, in Fig. HJwe assumed 
that r]/s = 0.2 is independent of -v/i. A recent hydrodynamic simulation suggests that rj/s 
is a decreasing function of ^/s [271], and this could alleviate the (small) discrepancy in the 
large-y^ region (remember that Av^ oc p/s). 


VII. SUMMARY AND CONCLUSIONS 


In this paper, we have revealed, in a completely analytical manner, a number of interesting 
features about the nature of hydrodynamics in the presence of conserved currents as well 
as the chemical potential (collision energy) dependence of the flow harmonics Let us 
summarize the main hndings. 


As a generalization of the Gubser flow 1^, we have derived an anisotropic solution 
of the relativistic Navier-Stokes equation coupled with conserved currents. Conformal 
symmetry and boost invariance have been assumed. The solution is valid to linear 
order in the shear viscosity rj and the eccentricity e„. Based on the solution and the 
Cooper-Frye formula, we analytically computed the flow harmonics Vn at hnite density. 


Multi-strange hadrons such as P and ^ may freeze out earlier than non-strange hadrons. Again this 
uncertainty mostly goes away in the ratio (15^ . 

It is worth mentioning that feed-down corrections (resonance decays) are not included in the current 


estimates and they could change the fitting parameters. 
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FIG. 4. Av^ = V2 

A) 

— V 

species of hadrons 

m, 

14| 


— V 2 {X) as a function of ^/s compared with the STAR data for five different 


color of data points is chosen to match the color of the corresponding line for each hadron species. 
We have used rj/s = 0.2, = 0.23//^ and /// = —0.15//^- 


• In ideal hydrodynamics, the QGP hreball follows a straight line trajectory /i/T = 
const, in the phase diagram in the (T, /i)-plane. The shear viscosity causes a deviation 
from the straight line as shown in fl32p . We expect this picture to be approximately 
correct in QCD in the deconhned phase of the hydro evolution. 


Vn is a decreasing function of density (or an increasing function of y/s) and decreases 
faster with n at higher densities Vn ~ g-’^in(4B3/27)_ 'jpjg jg because the lifetime of the 
hydrodynamic regime (~ ^/B^) is shorter at high density as it is correlated with the 
multiplicity ~ dN/dY through the constant energy condition (l20l) . In this regard, 


it is interesting to recall that in an early numerical study 


23|, a constant (or even 


decreasing) V 2 as a function of y/s was obtained if the hydro simulation is continued 
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to very low temperatures (the so-called ‘hydro limit’ The rising V 2 with en ergy 


24| 


can be obtained by switching off hydrodynamics at a relatively high temperature 
Our assumption of early freezeout is similar in spirit to this. 

At hnite chemical potential, there are new viscous corrections to Vn (the last terms 
proportional to 7 in fl^7|) ). Numerically, they are smaller than the contribution previ¬ 
ously found in the /i = 0 case |^. However, they give the leading order contribution 
to the difference in Vn between particles and antiparticles. 

The viscous corrections to are enhanced at high density. Even if rj/s is constant, 
the Knudsen number K grows at high density as it is inversely proportional to the 


multiplicity fl79p . At large-n, it is also enhanced linearly by n, Vn/v, 

Q. 


ideal 


1 - 0{nK) 


The elliptic flow difference between particles and antiparticles Avf = vf—V 2 (or more 
generally, Av^) can be nonzero only if the particle X is charged under some chemical 
potential(s) and the shear viscosity is nonvanishing. This is related to the deviation 
from constancy of the ratio /i/T due to viscous effects. Our result is summarized by 
the ‘master formula’ 0831) (or the more fundamental formula 0841 ) ) which schematically 
reads oc . This formula dictates the ordering Anf > Avf > Avf > Av^^ > 
0 > AnJ^ which seems to be borne out by the STAR result except for a few data 
points. Our mechanism of generating Avf is distinct from the previous theoretical 


considerations in 
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m 


44l | , but we hnd it has some common ground with the discussion 


431. Finally we pointed out that the observed magnitude of AnJ^ is large and can 


be htted only if we assume an unnaturally large value of the isospin chemical potential 
^j. This suggests that other mechanisms to generate Avif may be at work. 


Presumably some of the above features are empirically well known to the experts of hydro- 
dynamic simulations. However, they have not been systematically derived with the level of 
analytical detail presented in this work. 

There are a number of directions for future work. Admittedly, the assumptions of boost 
invariance and conformal invariance are too simplistic, especially at high density. One has 
to relax these approximations to be more realistic. Related to this, we only considered the 
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conformal equation of state e = 3p = T^f{fi/T) where the function / does not carry any 
information about the crossover and possibly hrst order phase transitions at hnite density. 
(Nevertheless it is remarkable that we can explain many features of Vn measured at different 
energies without such information.) It is important to hgure out how the presence of phase 
transitions in / is encoded in the observed behavior of Including the effects of anomaly 
(see, e.g., [4^) is also interesting. We hope to address these questions in future work. 
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Appendix A: Computation of (5J2,3 


In this Appendix, we carry out the computation of 572,3 dehned in fl54)) . We hrst note 
that the last two terms of fl5^ may be combined as 


/ , , X pt . , , ^ {‘^LfpT 

-Ptcos( 0 - 0p)--1-sm(0- 0p)— 


dx±_ x± 


dcf) B^{L‘^ + x‘if 


4xx cos 0 + 3e„ 


2Lx_l 

L? + 


-Ax±_ + n- 


L2 


x±_ 


' cos n((p + (pp) cos 0 + n - — sin n(0 + 0p) sin 0 

x± 


(Al) 


where we shifted the integration variable 0 as 0 ^ 0 + 0p. We then make the following 
replacement 

cos n(0 + 0p) cos 0 cos ncj) cos 0 cos ncpp —- cos(n — 1)0 cos n(f)p , 


sin u(0 + 0p) sin 0 —)■ sin U0 sin 0 cos ncpp - cos(?7, — 1)0 cos ncpp , (A2) 


where we neglected sinn0p which will vanish after the 0-integral, and also cos(?7,-|-1)0 which 
will lead to subleading terms In+i{z) ~ after the 0-integral compared to /„_i(z) ~ 
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Thus, flAlD effectively becomes 


{2LfpT 




4 x_l cos (j) + 3 e„ 


2Lx± 

L? + 


-2x±_ + 


a:_L 


cos(?7, — l)0cosn0j; 


.(A3) 


Similarly, we can write 5U defined in fl3T]) as 


(Su 

5u±_ cos 0 cos n(0 + 0p)- - sin 0 sin n(0 + 0^ 

x_i_ 


^ 1 
-)■ -- 


5m I 


2 Mj_o 

3z / 2 Lx_l 


5m0 

X_L 
n 


cos(n — l)0cosn0j; 


4a:_L V 


-2a;_L + 


nLr 

x±_ 


cos(n — 1)0 cos n(j). 


p ■ 


(A4) 


Using these simplifications, we get 


6J2 + SJ3 ~ e^e 


^^{2Lf'yKpT f ,, u 


+3 


2Lxj 


— 


LP' + x^ 

{2Lf-iKpT 

B^ 


d(dx±dcj) 


4 + ^ 


(L2 + x\y 


-2x±_ + 


nlP 


d(dx±d(j ) , n sn e 


X_L 

[/ I 2Lx± 


{L^ + xlf \L^ + x 


4:X±5U ( -^U + -^ + k] cos< 


COs(n — l)0COS?7,0p 




-2a;_L + 


x± 


3^ ( + ^ + cos 0 + 3 + fc 


cos(n — l)0cosn0p , (A5) 


where 


U 


rriT 

~T~ 

C 


cosh(C — Y) + z cos 0. 


(A6) 


Let us now define 


Yn{a) 


(iCd0e“^ cosn0 = A7rKo{amT/Tc)In{ctz) 


1 / (1 z \ 

47riLo(amT/T,)-(y) . 


(A7) 


























Using this we obtain 


w,+ sj, « f dx, 




(L^ + \L^ + X 


2Lx 


X 


■'1 


-2xx + 


nL/^ 


X_L 


y (|^ K - 2 + (^ + '=| y ..-2 


3 I UkL, + A^K„_, 


a=l 




Xx 
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3n 


2Lx± 

(L2 + 4)3 


-2a;_L + 


nL^ 

a;_L 


(n - 1)! V2 


n—1 


L 

4/ 


rrtr 


Ki{mT/T^) + {n- l)i^o("ir/7'c) + kKo{mT/Tc) 


, 3 n^(n- l)r(3n) 
3n-l r(4n) 


X 


mr 

[4/ 


Ki{mT/Tc) + (n - l)it'o("iT/Uc) + kKo{mT/Z) 


(AS) 


The correction to can be computed analogously to (|58l) . and the result is reported in 
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